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We have studied numerically the effect of quenched site dilution on a first order phase transition 
in three dimensions. We have simulated the site diluted three states Potts model studying in detail 
the second order region of its phase diagram. We have found that the v exponent is compatible with 
the one of the three dimensional diluted Ising model whereas the rj exponent is definitely different. 
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I. INTRODUCTION 

The study of the impurities effect on the critical be- 
havior of a pure material is an important issue since 
frequently real systems cannot be considered as pure. 
Nowadays the effect of dilution (disorder coupled to the 
energy density) on second order phase transitions is well 
understood. The Harris' criterionlj states that if the 
specific-heat of the pure system presents a power-like di- 
vergence (i.e. ttpuro > 0) the disorder induces a new Uni- 
versality Class. Otherwise (apuro < 0) the critical behav- 
ior of the model remains unchanged. The criterium does 
not decide in the marginal case ckpurc — 0. Moreover, it 
is possible to show rigorously that for all the continuous 
phase transitions in presence of disorder, the correlation- 
length critical exponent v ijerifies v > 2/d, d being the 
dimensionality of the spacea. 

When the pure model shows a first order phase tran- 
sition the situation is more complicated. However, there 
are a set of important results, both numerical and ana- 
lytical. For instance, Aizenman and Wehia showed rigor- 
ously that in two dimensions when introducing disorder, 
its conjugated density becomes a continuous function of 
the thermodynamic parameters. For instance, in a site 
diluted model the conjugate density is the energy, while 
in the Random Field Ising Model (RFIM) , it is the mag- 
netization. If the phase transition becomes continuous, 
one may ask about its Universality Class. A widely stud- 
ied model in this context has been the q-states Potts 
model, whose pure version in d = 2 undergoes a first 
order phase transition for q > 5. In recent numerical 
simulations using Monte Carlo or transfer matrix meth- 
ods, the v exponent has been found compatible within 
errors with the pure Ising value {v = 1) independently 
of q, but for the magnetic exponent, /?, the numerical 
results for flj pr 2 are significantly different from the pure 
Ising valueETEl. 

Unfortunately, for d > 2 the situation is not so clear. 



Cardy and JacobsenS (see also Cardyll3) have put forward 
a picture of the general behavior by means of a mapping 
from the diluted g-states Potts model to the RFIM. Their 
mapping being asymptotically exact for large q, their re- 
sults are also expected to hold for phase transitions with 
large latent heat. According to them, when a system 
that undergoes a first-order phase transition gets weakly 
diluted, the latent heat decreases. For larger dilutions, 
it will eventually vanish at a so-called tricritical point. 
Cardy and Jacobsen relate the latent heat with the order 
parameter of the RFIM, whose behavior is governed by 
the zero-temperature fixed point. In this way they are 
able to predict the critical exponents for the tricritical 
point P = /3rfim and v = i^rfim/(2 - aRFiM - /3RFiM)til- 
The critical-behavior of the system for larger dilutions 
remains unaddressed in their work. A physical realiza- 
tion of this scenario is provided by some magnetic semi- 
conductors like Zni„a;Mna;TeE3. The magnetic atoms of 
these materials behave as Heisenberg spins living in a 
fee lattice, with antiferromagnetic interactions. In these 
highly frustrated systems a first order phase transition 
is found in pure samples that gets second order upon 
dilution. 

We finish this overview describing the results obtained 
by Elderfield and Sherrington for the,diluted Potts Model 
in the Mean Field approximatiorO. They found that 
the phase transition is first order for q > 2, for all val- 
ues of spin-concentration. Another interesting and re- 
lated model is the Potts glass for g > 4 that, according 
to Mean Field theory, undergoes a first order transition 
with no latenl-heat while for q — 3 and 4 the transition 
is continuouscfl. 

In this work, we shall consider the effects of site- 
dilution in the three dimensional three state Potts model, 
whose pure version presents a weak first order transi- 
tion (small latent heat or large correlation length at the 
critical point). Our choice has been motivated by the 
ubiquity of weak first order transitions in nature. In par- 



ticular, the pure q — 3 Potts model shows very differ- 
ent experimental realizations appearing in very distant 
fields. For instance, we can cite the de-confining phase 
transition in quenched Quantum Chromodynamics. It 
can also characterize different situations in solid state 
physics. For instance a cubic ferromagnet with three 
easy axes of magnetization when a magnetic field in the 
diagonal of the cubic lattice is turned on (e.g. DyAl2), 
structural phase transitions (e.g. SrTiOa) ancLapme fluid 
mixtures of five (suitable chosen) componentalj. 

We now describe briefly the phase diagram of the three 
dimensional three states Potts model in the temperature- 
concentration plane {T,p) (see Fig. |l|). The pure model 
undergoes a (weak) first order phase transition, at a crit- 
ical temperature Tc{p — 1) separating the paramagnetic 
high-temperature phase from the low temperature or- 
dered one. This first order transition can be, in principle, 
continued inside the (T, p) plane, where the critical tem- 
perature Tc(p) will lower for smaller p. The latent-heat 
for the first-order phase transition will decrease until the 
tricritical point. At this point the model suffers a second 
order phase transition that continues (belonging to an- 
other Universality Class) until the Tc{pc) = percolation 
limit. We remark that, in the most economic picture, 
this phase diagram presents three different Universal- 
ity Classes: site percolation in three dimensions (which 
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has been studied in the literature, e.g. in Rcf. 
Universality Class of the tricritical point (conjectured in 
Ref. |l^) and the Universality Class that controls the crit- 
ical behavior in the region between the tricritical point 
and the percolation point. 

In this paper we will restrict ourselves to the study of 
the second-order line. Although an experimental real- 
ization of the site diluted Potts model is not yet known 
(disorder tends to couple with the order parameter rather 
than with the energy), whenever it will appear it will be 
interesting to have clear theoretical predictions at hand. 
The techniques used in this paper are well suited for sec- 
ond order transitions, but they should be modified in the 
concentration range for which the phase transition is first 
order. Work is in progress to study this regiontO. 



II. THE MODEL AND OBSERVABLES 

We have studied the three dimensional Site Diluted 
three state Potts Model, whose Hamiltonian defined on 
a cubic lattice with volume V — L^ \s 



n = Rc 
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and periodic boundary conditions are applied. In Eq. (|l|) 
ZiS are complex roots of z"^ = 1, and e^'s are uncorrelated 
quenched random variables which are 1 with probability 
p, and with probability 1 — p. The Boltzmann weight 
is proportional to exp(— /37i). 



We have used clusters algorithms in order to update 
the system. In a diluted system, the set of occupied sites 
can present regions that are lightly connected to the per- 
colating cluster. These regions are very dif&cjilt to equi- 
librate just with a single cluster algorithnJl3. We have 
found that a single cluster algorithm combined with a 
Heat Bath sweep per measure is efficient for large con- 
centrations. However, for small concentrations {p < 0.5) 
the previous method is not efficient enough due to the 
presence of intermediate-sized clusters, and we have used 
the Swendsen-Wang algorithmEj. 

We have simulated at p = 1.0, 0.9, 0.8, 0.7, 0.6, 0.5 and 
0.4 in lattices L = 8, 16, 32 and 64. For p = 0.8, 0.7, 0.4 
we have also run in i = 128 lattices. We have performed 
Nj = 200 nearly independent measures in every single 
disorder realization. For p < 0.8 the number of these 
realizations has been Ns = 10000, except for p = 0.8, 
L = 128, where we have fixed Ns = 1000. In the p = 0.9 
case we have measured in 2000 different disorder realiza- 
tions. The total amount of CPU time has been the equiv- 
alent of 16 years of 200 Mhz Pentium-Pro processor. For 
small diiutions we have performed the usual /3 extrap- 
olationEl while for p < 0.5 we used a p extrapolation 
methodEil. Let us recall that when planning a disordered 
model simulation, one should balance two competing ef- 
fects. First, to minimize statistical errors, it is better to 
work in a Nj ^ Ns regime. On the other hand, if Nj 
is too small, the usual calculation of f3 derivatives and 
extrapolations is biased. We follow the same procedure 
of Ref. ^ to eliminate the bias. With our simulation 
strategy {Nj <^ Ns), it is crucial to check that the sys- 
tem is sufficiently thermalized while taking measures. In 
order to ensure this, we have systematically compared 
the results coming from hot and cold starts: half of our 
statistics for the largest lattices have been obtained with 
hot starts while the other half comes from cold starts. 

Regarding the observables, in addition to the energy 
we have measured the complex magnetization 
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from which we obtain the real susceptibility as 

We have denoted with ((• • •)) the thermodynamical av- 
erage with fixed disorder and with (• ■ •) the average over 
the disorder. 

The formulae for the cumulants read 
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(74 being the standard Binder cumulant, (72 measures 
whether the susceptibility is or not a self-averaging quan- 
tity and (73 has been introduced since the three states 
Potts model is invariant under a global transformation 
of the Z3 group. The other cumulants, 32 and 174, are 
also trivially invariant since we have used the modulus of 
the complex magnetization in their construction. 

We have used a quotient methocO, in order to com- 
pute the critical exponents. We recall briefly the basis of 
this method. Let O be a quantity diverging as t~^° (t 
being the reduced temperature) in the thermodynamical 
limit. We can write the dependence of O on L and t in 
the following way 



0{L,t) =L^o/'' 



C(oo,t) 



+ o(i-,r 



, (7) 



where Fq is a (smooth) scaling function and (— w) is the 
biggest non positive eigenvalue of the Renormalization 
Group transformation (the corrections-to-scaling expo- 
nent). This expression contains the not directly measur- 
able term ^(00, t), but if we have a good definition of the 
correlation length in a finite box S,{L,t), Eq. (^ can be 
written 
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where Go is a smooth function related with Fq and F^ 
and we have neglected the term ^^ because we are work- 
ing deep in the scaling region. The definitioi^f the cor- 
relation length on a finite box that we use is 
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where x was defined in Eq. (j^) and F is given by 
V 
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z(fci, k2, fca) being the discrete Fourier transform of CiZi. 
We remark that the definition in Eq. (|^) makes sense as 
a correlation length only in the pure paramagnetic phase 
of the model. 

The main formula of the quotient method is 
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e.g. we compute the quotient between 0[sL,t) 
and 0{L,t) at the reduced temperature, t, in which 
(,{sL,t)/£^{L,t) — s. As particular cases of interest we 
cite the susceptibility, x, and the /3-derivative of the cor- 
relation length, df3^, whose associated exponents are: 



xa,,^ = 1 + 1^ , 
a^x = (2 - V)!^- 



(12) 
(13) 



A clean measure of scale invariance is provided by 
{^/L)\q ^^. Let us recall that ^/L is a monotonically 
growing function of the inverse temperature. In the or- 
dered phase it grows as L'^/^, while in the disordered 
phase decreases with growing lattice size. Therefore, for 
any pair of lattice sizes, there is a crossing temperature 
where Q^ = 2. In a second order transition, ^/L at the 
crossing point should tend to a non-vanishing universal 
value. For a first order transition, the crossing tempera- 
tures tend to the transition point but ^/ L at the crossing 
diverges due to the coexistence of ordered and disordered 
phases. 

We finally analyze the quotient of the cumulants 172, 
33 or 34 at two different lattices, L and sL, computed at 
the temperature where Q^ = s. Notice that for a second 
order phase transition the asymptotic limit (L —^ 00) 
of these quotients is 1 corrected by terms like L~'^ (see 
Eq. (|ll|)). ^ 

The quotient method, Eq. (pT|), has several interesting 
features. First, we profit of the large statistical corre- 
lation between Qo and Q^. Next, one does not need 
a previous estimate of the infinite volume critical point. 
Finally, it allows a simple control of the scaling correc- 
tions. All of this makes the method specially efficient for 
the measures of anomalous dimensions. 



III. NUMERICAL RESULTS 

The phase diagram that we have obtained numerically 
is shown in Fig. |lj and it agrees with the standard pic- 
ture. The dashed line corresponds to the part of the 
transition line where we do not find a clear second order 
asymptotic behavior. A clear first order signature is very 
difficult to see even for dilutions as small as p = 0.95. 
In this region the individual samples usually present a 
double peak structure. Consequently the thermalization 
process is very slow and a different algorithm for updat- 
ing muat be used. Work is in progress to analyze this 
regiontZl. 

Our scope is now to compute the critical exponents in 
the region in which the transition is clearly second or- 
der, i.e. the study of the Universality Class between the 
tricritical and percolation limits. The first stage is to de- 
termine where an asymptotic second order behavior has 
been reached with lattice sizes up to L = 128. In Fig. g 
we show the value of £,/L, at the points for which Q^ — 2 
for the different {L, 2L) lattice pairs and as functions of 

We have used fotA; the corresponding value of the site 
diluted Ising modeO. For p < 0.7 we find that ^/L seems 
to tend to a dilution-independent value. Notice the clear 
divergence for p = 1, where the transition is known to be 
first order. For p = 0.9 we find a similar trend that for 
the pure case while for p = 0.8 we find a transient be- 
havior: for small lattices ^/L grows while in the largest 
lattices it seems to approach the universal value. 
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FIG. 1. Phase diagram of the model defined in Eq. (|l|). 
The critical line has been drawn using the critical points ob- 
tained in our numerical simulations (the filled circles in the 
plot). 
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FIG. 3. The cumulant QaiL) at the points where 
5(2L,t) = 2^{L,t) as a function of L"""* for the different 
concentrations considered. 
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FIG. 2. The ratio ^/L as a function of L~ ' for different 
values of the concentration p. Notice that the behavior for 
p > 0.9 is clearly different from the rest of concentrations, 
where a non-vanishing, dilution-independent limiting value is 
likely to occur. 



Another interesting quantity in order to clarify the sec- 
ond order behavior is the cumulant g^. The g3{L) val- 
ues at the points where ^{2L,t) = 2^{L,t) are shown in 
Fig. 0. In this case we see a different scaling behavior 
for p = 0.9 and p — 0.8 up to the studied lattice sizes. 
We also guess from this figure that the uj value cannot 
be much larger than 0.4. 
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FIG. 4. Quotients of the cumulants 172, ga and g4 (filled 
squares, open squares and open circles, respectively) as a func- 
tion of L~ ' . Notice the different y scale in the p — 0.9 case. 

We have next considered the quotients of the different 
cumulants gi at the points where Q^ = 2. We recall that 
these quantities should go to 1 as i tends to infinity in a 
second order phase transition. We present our results in 
Fig. |j. At concentrations p — 0.9 and p = 0.8 we do not 
find an asymptotic behavior. For p = 0.7, the behavior 
is yet not monotonous. Only for p = 0.4,0.5 and 0.6 it 
seems that the asymptotic behavior is reached. Unfortu- 
nately, a reliable estimate of lo cannot be obtained but 
our results point to a value near 0.4. Moreover, the higher 
order scaling corrections are rather strong for these quan- 
tities. Finally, let us remark that the corrections to scal- 
ing and statistical errors are much larger for 52 and (74 
than for g^. 



Therefore, for the study the second order region, we 
conclude that only for p < 0.6 an asymptotic scaling be- 
havior for the considered lattice sizes has been found. 



0.9 



0.8 



0.7 



0.6 



0.5 



0.4 



8 0.571(4) 0.633(3) 0.662(4) 0.685(3) 0.706(4) 0.738(5) 

16 0.592(7) 0.659(3) 0.686(3) 0.692(3) 0.698(4) 0.711(5) 

32 0.664(12) 0.700(4) 0.695(4) 0.688(4) 0.694(4) 0.696(4) 

64 — 0.711(13) 0.707(4) — — 0.692(4) 

TABLE I. Apparent critical exponent v, obtained from 
Qdgi measured where Qj = 2 for all the concentrations stud- 
ied. 



L p^O.9 p^O.8 p^O.7 p = 0.6 p = 0.5 p = 0.4 

8 0.048(2) 0.057(2) 0.065(2) 0.0745(15) 0.079(4) 0.072(4) 

16 0.036(4) 0.045(2) 0.068(2) 0.0773(14) 0.079(2) 0.077(3) 

32 -0.029(9) 0.050(3) 0.074(3) 0.077(2) 0.079(2) 0.077(3) 

64 — 0.064(6) 0.071(3) — — 0.080(3) 

TABLE IL Apparent critical exponent rj, computed using 
Qy, for all the considered concentrations. 



We report the results for the critical exponents as func- 
tions of p and L in tables | and |l[ We have applied 
Eq. ( pdl ) with s = 2 to dpS, for computing z/ and to x 
for extracting rj. We can observe that the asymptotic 
behavior of these estimates for p > 0.7 is not clear. 
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FIG. 5. The v apparent critical exponent for three differ- 
ent concentrations (p = 0.6, 0.5 and 0.4) as a function of l/L. 
The shadowed region corresponds to-, the value of the three 
dimensional site diluted Ising modeEj. 
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FIG. 6. The finite size estimate for r\ critical exponent for 
three different concentrations (p = 0.6, 0.5 and 0.4) as a func- 
tion of 1/L. We show the value corpepponding to the three 
dimensional site diluted Ising modelESI with a shadowed re- 
gion. 

We have plotted the v and 77 apparent critical expo- 
nents for p — 0.6, 0.5 and 0.4 as functions of 1/L in 
Figs, a and P respectively. We recall that we have found a 
u: w 0.4 value for the cumulants. From the Figs, g and g 
we see that cu — 1 for the leading scaling-corrections term 
could be a reasonable choice in this case. A possible ex- 
planation of this contradiction could be that for the ob- 
servables used for computing the critical exponents the 
leading term (w w 0.4) vanishes. In any case we should 
remark that we have not a precise control over the scal- 
ing corrections unlike, for example, in the invesiigation of 
the three dimensional site diluted Ising modeEj. Fortu- 
nately, the scaling corrections for the critical exponents 
are rather small. Thus, it is not essential in this model 
to perform an infinite-volume extrapolation of our esti- 
mates. This is in marked contrast with the Ising case, 
where the extrapolation procedure was crucial to cor- 
rectly compute the critical exponents. 

From Figs. g| and ^ we estimate the critical exponents 
in the second order region as the displayed in the last 
row of the table |l|. We remark that the v value is in- 
distinguishable from the 3D site diluted Ising model one. 
However, the 77 value is very different from the values 
found for the rest of the models reported in the table II] . 



Model 



Pure IsingL 
Diluted Isiis 
PercolationtJ ■ — 1 
Tricritical PoinS 
This work 



0.6294(10) 


0.0374(12) 


0.3265(4) 


0.684(5) 


0.037(5) 


0.355(3) 


0.8765(18) 


-0.0460(3) 


0.4181(9) 


0.68(5) 


0.50(5) 


0.00(5) 


0.690(5) 


0.078(4) 


0.372(4) 



TABLE III. Critical exponents for some three dimensional 
models. 



IV. CONCLUSIONS 

We have numerically studied the three dimensional site 
diluted three states Potts model. The phase diagram in 
the temperature-concentration plane consists of a ferro- 
magnetically ordered phase separated from a paramag- 
netic, high temperature one. Between both regions there 
is a critical line, which is (weakly) first order in the limit 
of pure samples. For small concentrations, a clear second- 
order behavior is found, while the region with p ^ 0.9 
shows a different behavior, probably corresponding to a 
cross-over, more diflficult to analyze. 

The critical exponents of the second order region have 
been computed using Finite Size Scaling Techniques. We 
have found that these exponents are dilution indepen- 
dent, and that they show a very mild evolution with 
the lattice size. That is why a sound estimate of the 
critical exponents can be given, in spite of the fact that 
we have been unable of measuring the scaling-corrections 
exponent lu. This is in marked contrast with the situa- 
tion in the site-diluted Ising model, where the scaling- 
corrections are severe but uj can be obtained with a 15% 
accuracy. 

Regarding the variation of the critical exponents with 
q we have compared the results for the Potts model and 
the Ising case. We have found that the v exponent varies 
slowly (or perhaps remains unchanged) with the q value 
whereas the /3 or 77 exponents show a clear variation on 
this parameter. This picture strongly reminds the ob- 
tained in numerical simulations in two dimensions. 

The study of the first-order region and of the critical 
behavior in the neighborhood of the tricritical point re- 
quires rather different numerical techniques and will be 
the matter of future work. 
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